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Abstract 

Explicit formulas for the mean and variance of linear stochastic differential equa- 
tions are derived in terms of an exponential matrix. This result improved a pre- 
vious one by means of which the mean and variance are expressed in terms of 
a linear combination of higher dimensional exponential matrices. The important 
role of the new formulas for the system identification as well as numerical algo- 
rithms for their practical implementation are pointed out. 
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1. Introduction 

In previous works [8, 9], explicit formulas for the mean and variance of linear 
stochastic differential equations (SDEs) where derived in terms of a linear com- 
bination of exponential matrices. The formulas were obtained as solution of the 
ordinary differential equations (ODEs) that describe the evolution of the mean and 
variance. Apart from being of mathematical interest for control engineers, these 
formulas have played a crucial role in the practical implementation of suboptimal 
linear filters [8], Local Linearization filters [9] and approximate Innovation esti- 
mators [10] for the identification of continuous-discrete state space models. In 
a variety of applications, these methods have shown high effectiveness and effi- 
ciency for the estimation of unobserved components and unknown parameters of 
SDEs given a set of discrete observations. Remarkable is the identification, from 
actual data, of neurophysiological, financial and molecular models among others 
(see, e.g., [1, 3, 5, 13, 14]). Therefore, a simplification of the formulas for the 
mean and variance of linear SDEs might imply a sensible reduction of the compu- 
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tational cost of the mentioned system identification methods and, consequently, a 
positive impact in applications. 

In this paper, simplified explicit formulas for the mean and variance of linear 
SDEs are obtained in terms of just one exponential matrix of lower dimension- 
ality. The formulas are derived from the solution of the ODEs that describe the 
evolution of the mean and the second moment of the SDEs. The variance is then 
obtained from the well-known formula that involves the first two moments. The 
computational benefits of the simplified formulas are pointed out. 

2. Notation and Preliminaries 

Let us consider the J-dimensional linear stochastic differential equation 

m 

dx(t) = (Ax(t) +a{t))dt + £(B ; x(?) +b i {t))dw i {t) (1) 

for all t e [to, T] , where w = (w 1 , . . . , w m ) is an m-dimensional standard Wiener 
process, A and B, are d x d matrices, and a(f) = ao + &it and b,-(f) = b^o + b^i* 
are J-dimensional vectors. Suppose that there exist the first two moments of x for 
all t G [to,T\. 

The ordinary differential equations for the J-dimensional vector mean m t — 
E(x(t)) and the d x d matrix second moment P t = E(x(t)x J (t)) of x(t) are [11] 

dm, 

— =Anw + .(0 

and 

dP m 
— - = AP f + P f A T + £ B^BJ +#(*), 
at (=1 

where 

m 

®{t) =a(Om f T + m,a T (0 + £B i ni f b ; T (0+b,-niTBT(0+b i (ObJ(0. (2) 

The solution of these equations can be written as [8, 6] 

m, = m + Le c( ^ o) r (3) 



2 



and 



f-fo 



vec(P t ) = e-^-'o) (vec(Po) + j e-^ s vec(^(s + t ))ds) , (4) 

o 

where mo = E(x(to)) and Po = E(x(to)x J (to)) are the first two moments of x at 
to, and the matrices C, L and r are defined as 



C = 



A ai Am + a(? ) 
1 




G <%(d+2)x(d+2) 



(5) 



L = [ Lj 0^x2 ] and r T = [ lx ( rf+1 ) 1 ] for non-autonomous equations (i.e., 
with non zero ai, b ; .i); and as 



A Am + a(? ) 




gj(rf+l)x(d+l)_ 



(6) 



L= [ Lj Odxi ] andr T = [ Oi X d 1 ] for autonomous equations (i.e., with ai = 
b; i = 0). Here, 



^/ = a©a + £b ; -®bJ 

(=i 



(7) 



is a J 2 x d 2 matrix, and Id is the J-dimensional identity matrix. The symbols 
vec, © and © denote the vectorization operator, the Kronecker sum and product, 
respectively. 

The following lemma provides simple expressions to computing multiple in- 
tegrals involving matrix exponentials such those appearing in (4). 

Lemma 1. ([15]) Let Ai, A 2 , A3 and A4 be square matrices, n\, n 2 , n$ andn^ be 
positive integers, and set m to be their sum. If the mxm block triangular matrix 
M is defined by 



M 



Ai 


Bi 


Ci 


Di " 


} n\ 





A 2 


B 2 


c 2 


} n 2 








A3 


B 3 


} m 











A 4 


} «4 



3 



then for s >0 



Fi(j) Gi(s) Hi(j) Ki(s) 

F 2 (s) G 2 (j) H 2 (s) 

F 3 (s) G 3 (s) 

F 4 (s) 



exp(sM), 



where 

Fj(s)=exp(Ajs),forj= 1,2,3,4 

s 

Gj(s) = J exp(Aj(s — u))Bjexp(Aj + iu)du,for j= 1,2,3 


s 

Hj(s) = J exp(Aj(s — u))Cjexp(Aj + 2u)du 
o 

+ J Jexp(Aj(s — u))Bjexp(Aj + i(u — r))Bj + iexp(Aj + 2r)drdu 



.v 

Ki(*) = Jexp(Ai(s — u))Diexp(A4ii)du 


+ exp(A 1 (^ — M))[Ciexp(A 3 (M — r))B 3 + B 1 exp(A 2 (i/ — r))C2]exp(A 4 r)Jr<ii<. 
o o 

.? m r 

J J J exp(Ai(s — u))Biexp(A2(u — r))B2exp(Ai,(r — w))BT l exp(A4w)dwdrdu. 



+ 

boo 

A generalization of the above lemma for integrals with higher multiplicity is 
given by Theorem 1 in [2], which will be applied as well in a proof below. 

3. Simplified formulas for the first two moments 

In this section simplified formulas for the first two moments of the linear SDE 
(1) and two of their special forms are derived. Equations with additive and mul- 
tiplicative noise as well autonomous and nonautonomous equations are distin- 
guished. 
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3.1. Equations with multiplicative noise 
Lemma 2. 

vec(SS(s + t )) = ^i+ M 2 s + 3§ 3 s 2 + ^ 4 e Ci r + s& 5 e Cs r, (8) 

for all s>0, where the vectors SS\ = vec{fi\) + famo, S$2 = veci fa) + P5 m o and 
SB?, = vec(fi?,), and the matrices = J84L and = J65L are defined in terms of 
the matrices 

m 

Pi = X)(bi,o+bi,ifo)(bi,o+bi,if ) T 

!=1 

m 

Pi = £ (b,,o + b„ i*o)bj! + b,-, 1 (b ; - + b M * ) T 

m 

J3 4 = (a + ai?o)©(ao + ai?o) + J2(bi )0 +b ( - ) i?o) ®B;+B i <8>(b i> o + b /i ifo) 

1=1 

m 

jS 5 = ai©ai + £bi > i®B i +B i ®b l - > i. 

iVoo/ The formula for vec(SS(s + to)) is directly obtained by substituting (3) in 
(2) witht = s + t . □ 

The main result of this paper is the following. 

Theorem 3. Let x be the solution of the linear SDE (1) with moments mo = 
E(x(to)) and Pq = E (x(to)x T (to)) at to. Then, the first two moments ofx can 
be computed as 

m f =m + L 2e M ^u 

and 

vec(P t )=L l e M{t - t ^u 
for all t G [to, T], where the (d 2 + 2d + 7) —dimensional vector u and the matrices 
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M, Li, L2 are defined as 





" srf 


^5 


^4 




^2 


m x ' 




vec(Po) 







c 


















M = 








c 











, u = 


r 














2 

























1 





























1 



and 

L 2 = [ ®dx{d 2 +d+2) Id ®dx5 ] , Li = [ l d 2 d 2 x ( 2d+7 ) ] , 

with £$[ defined as in Lemma 2, C, r in (5), and srf in (7). 
Proof. From (4) and (8) follows that 

vec(P t ) = F!vec(P ) + Ki + H ir 

where 



K, 



and 



j e ^(t-to-s)gg lds+ J J e ^^- tQ - s ^ 2 duds + 2 



t-t 



SS^drduds 







000 



H, 



t— to t—t(, s 

J e ^ t - t °- s ^ 4 e Cs ds+ j je^-^SSse^duds. 







Further, with F3 = e c ( r to \ (3) can be written as 

m t = m + LF 3 r, 

where the matrices L, C, r are defined as in (5). Thus, by a direct application 
of Theorem 1 in [2] (alternatively, Lemma 1 for Fi, F3 and Hi can be used by 
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simplicity) follows that 



Fx - Hi - - 

- - - - 

F 3 - - 

- - 

- 





Ki 




e M ^ withM: 



at 








^5 ^4 

c 





















2 










1 





where 38[ are defined as in Lemma 2, and st in (7). This implies that 

m, = mo + Lv2 

and 



vec(P t ) 



vi, 



where the J -dimensional vector vi and the J-dimensional vector V2 are defined 

as 

vec(Po) 

r 


1 



v 2 



e M (t-to) u with u 



Proof concludes by verifying that Lv2 = L2V and Vi = Liv, where v T = [vi, — , V2, 

-],L 2 = [ dx(d 2 +rf+2 ) L dx3 ] andLi = [ l d i ^ d ^x{2d+i) ]■ □ 

For autonomous equations with multiplicative noise the formulas of the previ- 
ous theorem can be simplified as follows. 

Theorem 4. Let x be the solution of the linear SDE (1) with ai = b; 1 = 0. Let 

mo = E(x(to)) and Pq = E (x(to)x J (to)) be moments ofx at to. Then, the first two 
moments ofx can be computed as 



m t = m + L 2 e 



M(f-fo) 



u 



and 



vec(P t )=Ue M{t - to) u 



for all t e [to, T], where the (d 2 + d + 2)— dimensional vector u and the matrices 
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Mi, Li are defined as 





" sf 


9S\ 


^4 " 




vec(Po) 


M = 











, u = 


1 










c 




r 



and 

L 2 = [ °dx(d 2 +l) l d Orfxl ] , Ll = [ l d 2 d 2 x ( d+2 ) ] 

with 3§\ and S3 a, defined as in Lemma 2, C, r in (6), and srf in (7). 

^ 5 = 0. Thus, from (4) and (8) follows 

vec(P t ) = F!vec(P ) +Gi +H ir , 



Proof. Since ai = b f) i = 0, ^ 2 = ^3 
that 



where 



G, 



t-t 



SS\ds 



and 



o 

t-to 



H, 



o 



By a direct application of Lemma 1 follows that 

m f = mo + LF3r 

and 

vec(P t ) = Fivec(Po) +Gi +Hir, 
where the matrices Fi, F3, K| , and Hi are defined as 



Fi Gi Hi 
- - 
F 3 



e M(t-t ) with M 




c 



with SS\ and 884 defined as in Lemma 2, C, r in (6), and s& in (7). This implies 
that 

m, = mo + LV2 
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and 



vec(P t ) = vi. 



where the J -dimensional vector vi and the d-dimensional vector V2 are defined 

as 



vi 

v 2 



e M(r-f ) u withu 



vec(Po) 
1 
r 



Proof concludes by verifying that Lv 2 = L 2 v and vi = Liv, where v T = [vi, — , V2], 
L2 = [ dx(d 2 +1) L ] andLi = [ l d i rf 2 x(d+2) ]. □ 

3.2. Equations with additive noise 

For autonomous SDEs with additive noise an additional simplification of the 
explicit formulas for the first two moments can be archived. 

Theorem 5. Let x be the solution of the linear SDE (1) with B, = and ai = 
b J; i = 0. Let mo = E(x(to)) and Po = E(x(to)x T (to)) be moments ofx at to. Then, 
the first two moments ofx can be computed as 



m, = mo + ki 



and 



P, = FiPoFj+HiFj+FiHj 

for all t e [to, T], where the matrices F\, Hi, and ki are defined as 

Hi ki 



Fi 


0- 




Mt-to) 



being 



M 



A 




m 

a mj + i Eb f , bT 

j=i 

(Am + a ) T 


Am + a 

















-AT 


















a (2J + 2) x (2J + 2) matrix. 
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Proof. Thanks to B ; = 0, ai = bn = 0, the commutation of the matrices (I® A) 
and (A® I), and by using expressions (4) and (2) it is obtained that 

t-to 

p t = e A(t-to)p oe AT(t-t ) + J e A( f - ?0 -^( s + ?0 ) e AT( f - ?0 -,)^ 



where now 38 {s + ?o) reduces to 

m 

m(s + 1 ) = a mj +f0 + m s+to al + £ b,-,obT . 

i=\ 

Taking into account that m s+to can be rewritten as [8] 

m s+tQ =m + ki, 

with 

.? 

k! = Je A{s - ll) (Amo + a )du 
o 

we have that 

s 

a mj +f() = a mj +ao(Am + ao) T y e AJ( - s ^du, 

o 

and so P f can be rewritten as 

P f = F 1 P F{+H 1 F{ + F 1 H{, 

where 

Fi = e A ('-<^ 
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and 



Hi 



/ 



e Mt-to 



i m 

5) (aom5 + -£b i)0 bJ ) e - AT ^ 



o 



1=1 



.v 



+ 



/ 



^aolAmo + ao) 1 e AJ "duds. 



o 



o 



Proof concludes by a direct application of Lemma 1 



□ 



4. Computational aspects 

Theorem 3 provides explicit formulas for the first two moments of the linear 
SDE (1) in terms of an exponential matrix of dimension d 2 + 2d + l. By using the 
well-known expression 



for the variance var(z) of a random variable z, the variance var(x) of x solution 
of (1) can be straightforward computed as 



for all t e [to, T], where m t and P ? are given as in Theorem 3. 

By taking into account that the explicit formulas for the mean and variance 
of (1) obtained in [8, 9] involve the computation of seven exponential matrices of 
different dimensions up to a maximum of 3d 2 + Ad + 4, it is obvious the remark- 
able advantages of the new simplified formulas. From a computational viewpoint, 
this includes a considerable reduction of the computer storage capacity and the 
computational time required for their evaluations through the well known Pade 
method [12, 7] for exponential matrices. But, in addition, the new formulas allow 
the efficient computation of the mean and variance of high dimensional systems 
of the linear SDEs by means of the Krylov subspace method [12, 7] for exponen- 
tial matrices, which is crucial in many practical situations. Other advantage of the 
exponential form of these formulas is the flow property the exponential operator, 
which allows an extra reduction of the computational time when the mean and 
variance of (1) are required on consecutive time instants with multiplicity. In this 
case, the firsts two conditional moments at the first time instant after the initial 



var(z) = £(zz T ) -£(z)£(z T ) 



var(x(t)) = P t -m t m] 
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one is computed through the exponential matrix of Theorem 3, whereas the others 
at the remainder times are computed by simple multiplications of the exponential 
matrix just mentioned. 

Theorems 4 and 5 provide explicit formulas for the first two moments of au- 
tonomous linear SDEs, which involve an exponential matrix of lower dimension- 
ality: d 2 + d + 2 for equation with multiplicative noise, and 2d + 2 for equations 
with additive noise. This yields extra advantages in a number of important appli- 
cations. 

The formulas derived in this paper are one of the main components of the 
adaptive Local Linearization filters [4], which provide accurate and computation- 
ally efficient approximations to the solution of the linear minimum variance filter- 
ing problem for continuous-discrete state space models. 
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